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Abstract 

Calculating the structure of white dwarf and neutron stars would be a suitable topic for an 
undergraduate thesis or an advanced special topics or independent study course. The subject is 
rich in many different areas of physics accessible to a junior or senior physics major, ranging from 
thermodynamics to quantum statistics to nuclear physics to special and general relativity. The 
computations for solving the coupled structure differential equations (both Newtonian and general 
relativistic) can be done using a symbolic computational package, such as Mathematica. In doing 
so, the student will develop computational skills and learn how to deal with dimensions. Along the 
way he or she will also have learned some of the physics of equations of state and of degenerate 
stars. 

PACS numbers: 01.40.-d, 26.60.+C, 97.10.Cv, 97.20.Rp, 97.60.Jd 
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I. INTRODUCTION 



In 1967 Jocelyn Bell, a graduate student, along with her thesis advisor, Anthony Hewish, 
discovered the first pulsar, something from outer space that emits very regular pulses of 
radio energy. After recognizing that these pulse trains were so unvarying that they would not 
support an origin from LGM's (Little Green Men), it soon became generally accepted that 
the pulsar was due to radio emission from a rapidly rotating neutron star 1] endowed with a 
very strong magnetic field. By now more than a thousand pulsars have been catalogued ^ . 
Pulsars are by themselves quite interesting Q], but perhaps more so is the structure of the 
underlying neutron star. This paper discusses a student project dealing with that structure. 

While still at MIT before coming to Los Alamos, one of us (Reddy) had the pleasure of 
acting as mentor for a bright British high school student, Aiden J. Parker. Ms. Parker was 
spending the summer of 2002 at MIT as a participant in a special research program (RSI). 
With minimal guidance she was able to write a Fortran program for solving the Tolman- 
Oppenheimer-Volkov (TOV) equations [J] to calculate masses and radii of neutron stars 
(!)• 

In discussing this impressive performance after Reddy' s arrival at LANL, the question 
came up of whether it would have been possible (and easier) for her to have done the com- 
putation using Mathematica (or some other symbolic and numerical manipulation package). 
This was taken up as a challenge by the other of us, who also figured it would be a good 
opportunity to learn how these kinds of stellar structure calculations are actually done. (Sil- 
bar's only previous experience in this field of physics consisted of having read, with some 
care, the chapter on stellar equilibrium and collapse in Weinberg's treatise on gravitation 
and cosmology 0.) 

In the process of meeting the challenge, it became clear to us that this subject would 
be an excellent topic for a junior or senior physics major's project or thesis. After all, if a 
British high school student could do it. . . There is much more physics in the problem than 
just simply integrating a pair of coupled non-linear differential equations. In addition to 
the physics (and even some astronomy), the student must think about the sizes of things he 
or she is calculating, that is, believing and understanding the answers one gets. Another 
side benefit is that the student learns about the stability of numerical solutions and how 
to deal with singularities. In the process he or she also learns the inner mechanics of the 
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calculational package (e.g., Mathematica) being used. 

The student should begin with a derivation of the (Newtonian) coupled equations, and, 
presumably, be spoon-fed the general relativistic (GR) corrections. Before trying to solve 
these equations, one needs to work out the relation between the energy density and pressure 
of the matter that constitutes the stellar interior, i.e., an equation of state (EoS). The first 
EoS's to try can be derived from the non-interacting Fermi gas model, which brings in 
quantum statistics (the Pauli exclusion principal) and special relativity. It is necessary to 
keep careful track of dimensions, and converting to dimensionless quantities is helpful in 
working these EoS's out. 

As a warm-up problem the student can, at this point, integrate the Newtonian equations 
and learn about white dwarf stars. Putting in the GR corrections, one can then proceed in 
the same way to work out the structure of pure neutron stars (i.e., reproducing the results 
of Oppenheimer and Volkov ^). It is interesting at this point to compare and see how 
important the GR corrections are, i.e., how different a neutron star is from that which 
would be given by classical Newtonian mechanics. 

Realistic neutron stars, of course, also contain some protons and electrons. As a first 
approximation one can treat this mult i- component system within the non-interacting Fermi 
gas model. In the process one learns about chemical potentials. To improve upon this 
treatment we must include nuclear interactions in addition to the degeneracy pressure from 
the Pauli exclusion principle that is used in the Fermi gas model. The nucleon-nucleon 
interaction is not something we would expect an undergraduate to tackle, but there is 



a ..p,e .ode, ,wMc. we ,ea™d a.o. .o. P.a^.h ft fo. ..e „.c,ea. .a„e. EoS. 

It has parameters which are fit to quantities such as the binding energy per nucleon in 
symmetric nuclear matter, the so-called nuclear symmetry energy (it is really an asymmetry) 
and the (not so well known) nuclear compressibility. Working this out is also an excellent 
exercise, which even touches on the speed of sound (in nuclear matter). With these nuclear 
interactions in addition to the Fermi gas energy in the EoS, one finds (pure) neutron star 
masses and radii which are quite different from those using the Fermi gas EoS. 

The above three paragraphs provide the outline of what follows in this paper. In the 
presentation we will also indicate possible "gotcha's" that the student might encounter and 
possible side-trips that might be taken. Of course, the project we outline here can (and 



probably should) be augmented by the faculty mentor 



with suggestions for byways that 
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might lead to publishable results, if that is desired. 

We should point out that there is a similar discussion of this subject matter in this journal 



by Balian and Blaizot |^|. They, however, used this material (and other, related materials) 
to form the basis for a full-year course they taught at the Ecole Polytechnique in France. 
Our emphasis is, in contrast, more toward nudging the student into a research frame of mind 
involving numerical calculation. We also note that much of the material we discuss here is 
covered textbook bv Shapi.o and TeuUols.y g Howeve. a. the .eade. will .ot.ee. 

the emphasis here is on students learning through computation. One of our intentions is 
to establish here a framework for the student to interact with his or her own computer 
program, and in the process learn about the physical scales involved in the structure of 
compact degenerate stars. 

II. THE TOLMAN-OPPENHEIMER-VOLKOV EQUATION 
A. Newtonian Formulation 

A nice first exercise for the student is to derive the following structure equations from 

classical Newtonian mechanics, 

dp _ Gp{r)M{r) _ Ge{r)M{r) 
dr c^r^ 

^ = 4.rV(r) = ^^Ai ^ (2) 

Mir) = 4n r r'^dr'ph') = An f r'^dr' eir') / . (3) 
JO Jo 

Here G = 6.673 x 10~^ dyne-cm^/g^ is Newton's gravitational constant, p{r) is the mass den- 
sity at the distance r (in gm/cm^), and e is the corresponding energy density (in ergs/cm^) 
|lO|. The quantity A4{r) is the total mass inside the sphere of radius r. A sufficient hint for 
the derivation is shown in Fig. ^ (Challenge question: the above equations actually hold 
for any value of r, not just the large-r situation depicted in the figure. Can the student also 
do the derivation in spherical coordinates where the box becomes a cut-off wedge?) 

Note that, in the second halves of these equations, we have departed slightly from New- 
tonian physics, defining the energy density e(r) in terms of the mass density p(r) according 
to the (almost) famous Einstein equation from special relativity, 

e(r) = p(r)c^ . (4) 



p(r+dr) = F(r+dr)/A 
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FIG. 1: Diagram for derivation of Eq. 

This allows Eq. ^ to be used when one takes into account contributions of the interaction 
energy between the particles making up the star. 

In what follows, we may inadvertently set c = 1 so that p and e become indistinguishable. 
We'll try not to do that here so students following the equations in this presentation can 
keeping checking dimensions as they proceed. However, they might as well get used to this 
often-used physicists' trick of setting c = 1 (as well as = 1). 

To solve this set of equations for p(r) and hAir) one can integrate outwards from the 
origin (r = 0) to the point r = R where the pressure goes to zero. This point defines R as 
the radius of the star. One needs an initial value of the pressure at r = 0, call it po, to do 
this, and R and the total mass of the star, Ai{R) = M, will depend on the value of po- 

Of course, to be able to perform the integration, one also needs to know the energy 
density e(r) in terms of the pressure p{r). This relationship is the equation of state (EoS) 
for the matter making up the star. Thus, a lot of the student's effort in this project will 
necessarily be directed to developing an appropriate EoS for the problem at hand. 
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B. General Relativistic Corrections 

The Newtonian formulation presented above works well in regimes where the mass of the 
star is not so large that it significantly "warps" space-time. That is, integrating Eqs. (1) 
and (2) will work well in cases when general relativistic (GR) effects are not important, such 
as for the compact stars known as white dwarfs. By creating a quantity using G that has 
dimensions of length, the student can see when it becomes important to include GR effects. 
(This happens when GM/c^R becomes non-negligible.) As the student will see, this is the 
case for typical neutron stars. 

It is probably not to be expected that an undergraduate physics major derive the GR 
corrections to the above equations. For that, one can look at various textbook derivations 
of the Tolman-Oppenheimer-Volkov (TOV) equation . It should suffice to simply state 
the corrections to Eq. (Q) in terms of three additional (dimensionless) factors, 

dp Ge{r)M{r) 
dr c^r"^ 

The differential equation for Ai{r) remains unchanged. The first two factors in square 
brackets represent special relativity corrections of order v'^/c^. This can be seen in that 
pressure p goes, in the non-relativistic limit, like kp/2m = mv'^/2 (see Eq. ()14|) below) while 
e and Aic^ go like mc^. That is, these factors reduce to 1 in the non-relativistic limit. (The 
student should have, by now, realized that p and e have the same dimensions.) The last 
factor is a GR correction and the size of GAi/c^r, as we emphasized above, determines 
whether it is important (or not). 

Note that the correction factors are all positive definite. It is as if Newtonian gravity 
becomes stronger at every r. That is, special and general relativity strengthens the relentless 
pull of gravity! 

The coupled non-linear equations for p{r) and Ai{r) can also in this case be integrated 
from r = for a starting value of Pq to the point where p{R) = 0, to determine the star mass 
M = Ai{R) and radius R for this value of po. These equations invoke a balance between 
gravitational forces and the internal pressure. The pressure is a function of the EoS, and for 
certain conditions it may not be sufficient to withstand the gravitational attraction. Thus 
the structure equations imply there is a maximum mass that a star can have. 



1 + 
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A^(r)c2 



2GM{r) 
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— 1 



(5) 
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III. WHITE DWARF STARS 



A. A Few Facts 

Let us violate (in words only) the second law of thermodynamics by warming up on 
cold compact stars called white dwarfs. For these stellar objects, it suffices to solve the 
Newtonian structure equations, Eqs. (dJ-Q [ill- 

White dwarf stars [l^ were first observed in 1844 by Friedrich Bessel (the same person 
who invented the special functions bearing that name). He noticed that the bright star 
Sirius wobbled back and forth and then deduced that the visible star was being orbited by 
some unseen object, i.e., it is a binary system. The object itself was resolved optically some 
20 years later and thus earned the name of "white dwarf." Since then numerous other white 
(and the smaller brown) dwarf stars have been observed (or detected). 

A white dwarf star is a low- or medium-mass star near the end of its lifetime, having 
burned up, through nuclear processes, most of its hydrogen and helium forming carbon, 
silicon or (perhaps) iron. They typically have a mass less than 1.4 times that of our Sun, 
Mq = 1.989 X 10^^ g They are also much smaller than our Sun, with radii of the order 
of 10^ km (to be compared with Rq = 6.96 x 10^ km). These values can be worked out from 
the period of the wobble for the dwarf-normal star binary in the usual Keplerian way. As a 
result (and as is also the case for neutron stars), the natural dimensions for discussing white 
dwarfs are for masses to be in units of solar mass, M©, and distances to be in km. Using 
these numbers the student should be able to make a quick estimate of the (average) densities 
of our Sun and of a white dwarf, to get a feel for the numbers that he will be encountering. 

Since GM/c^R ^ 10~^ for such a typical white dwarf, we can concentrate here on solving 
the non-relativistic structure equations of Sec. 2.1. (Question: why is it a good approxima- 
tion to drop the special relativistic corrections for these dwarfs?) 

The reason a dwarf star is small is because, having burned up all the nuclear fuel it 
can, there is no longer enough thermal pressure to prevent its gravity from crushing it 
down. As the density increases, the electrons in the atoms are pushed closer together, which 
then tend to fall into the lowest energy levels available to them. (The star begins to get 
colder.) Eventually the Pauli principle takes over, and the electron degeneracy pressure 
(to be discussed next) provides the means for stabilizing the star against its gravitational 
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attraction [9|, |13|. This is the physics behind the EoS which one needs to integrate the 
Newtonian structure equations above, Eqs. (jT]) and (j21). 



B. Fermi Gas Model for Electrons 

For free electrons the number of states dn available at momentum k per unit volume is 

(This is a result from their modern physics course that students should review if they don't 
remember it.) Integrating, one gets the electron number density 



n = / k'dk = . (7) 



{27ihy Jo Stt^^ 
The additional factor of two comes in because there are two spin states for each electron 
energy level. Here kpc, the Fermi energy, is the maximum energy electrons can have in the 
star under consideration. It is a parameter which varies according to the star's total mass 
and its history, but which the student is free to set in the calculations he or she is about to 
make. 

Each electron is neutralized by a proton, which in turn is accompanied in its atomic 
nucleus by a neutron (or perhaps a few more, as in the case of a nucleus like ^^Fe26)- Thus, 
neglecting the electron mass rrie with respect to the nucleon mass rriN, the mass density of 
the star is essentially given by 

p = nruNA/Z , (8) 

where A/Z is the number of nucleons per electron. For ^^C, A/Z = 2.00, while for ^^Fe, 
A/Z = 2.15 . Note that, since n is a function of kp, so is p. Conversely, given a value of p, 

\mN AJ 

The energy density of this star is also dominated by the nucleon masses, i.e., e ~ pc^. 

The contribution to the energy density from the electrons (including their rest masses) is 



6elec(M = J^r^^'^'+^yy^'^'dk 
/•kp/nieC 

= eo / {u^ + l)^/VrfM 
Jo 



(2x^ + x)(l + x^)^/2-sinh-^(x)l , (10) 



where 



4 5 



(11) 



carries the desired dimensions of energy per volume and x = kp/rrieC. The total energy 
density is then 

e = nmNA/Z + ecicc{kF) ■ (12) 

One should check that the first term here is much larger than the second. 

To get our desired EoS, we need an expression for the pressure. The following presents a 
problem (!) that the student should work through. From the first law of thermodynamics, 
dU = dQ — pdV, then at temperature T fixed at T = (where dQ = since dT = 0) 



p = - 



dU_ 

dV 



T=0 



^die n) de 

n — ; = n- e = na — e 

dn dn 



(13) 



where the energy density here is the total one given by Eq. (fT^ . The quantity /ij = de/drii 
defined in the last equality is known as the chemical potential of the electrons. This is a 
concept which will be especially useful in Section 5 where we consider an equilibrium mix 
of neutrons, protons and electrons. 

Utilizing Eq. (fTUI). Eq. (fT^ yields the pressure (another problem!) 



Stt 



t\k^c^ + mld')-^'^k^dk 
Jo 



fo 
3 

24 



kp/m^c 



+ l)-^/\^du 



{2x^ - 3x)(l + x^/^ + 3sinh~^l 



x\ 



(14) 



(Hint: use the n'^d{e/n) /dn form and remember to integrate by parts.) 

Using Mathematica [14j the student can show that the constant eo = 1.42 x 10^^ in units 
that, at this point, are erg/cm'^. (Yet another problem: verify that the units of eo are as 
claimed |l5|.) One also finds that Mathematica can perform the integrals analytically. (We 
quoted the results already in the equations above.) They are a bit messy, however, as they 
both involve an inverse hyperbolic sine function, and thus are not terribly enlightening. It 
is useful, however, for the student to make a plot of e versus p (such as shown in Fig. ^ for 
values of the parameter < kp < 2me. This curve has a shape much like e^^^ (the student 
should compare with this), and there is a good reason for that. 
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FIG. 2: Relation between pressure p (y-axis) and energy density e (x-axis) in the free electron 
Fermi gas model. Units are ergs/cm'^. Note that the pressure is much smaller than the energy 
density, since the latter is dominated by the massive nucleons. 

Consider the (relativistic) case when kp S> mg. Then Eq. fjl4|) simphfies to 

p{Kp) = — / u au = — [kp rrieC) = 

^ 3 Jo 12^ ' ' 127r2 \^ tunA ) 



(15) 



where 



he 



4/3 



rel 



127r2 {ArriNC^ J ' ^^^^ 
A star having simple EoS like p = Ke"' is called a "polytrope" , and we therefore see that the 
relativistic electron Fermi gas gives a polytropic EoS with 7 = 4/3. As will be seen in the 
next subsection, a polytropic EoS allows one to solve the structure equations (numerically) 
in a relatively straight-forward way lfi| . 

There is another polytropic EoS for the non-interacting electron Fermi gas model corre- 
sponding to the non-relativistic limit, where kp <^ rrie- In a way similar to the derivation of 
Eq. (fT3j). one finds 



P 



non— re. 



5/3 



where K, 



37r2Z 



5/3 



non— rcl 



(17) 



[Question: what are the units of i^Trei and -ft'non-rei? Task: confirm that in the appropriate 
limits, Eqs. (fTUj) and reduce to those in Eqs. (fT3j) and (|T7j) .] 
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C. The Structure Equations for a Polytrope 



As mentioned earlier, we want to express our results in units of km and Mq. Thus it is 
useful to define Ai{r) = A^(r)/M0. The first Newtonian structure equation, Eq. (Q), then 
becomes 

where the constant Rq = GMq/c^ = 1.47 km. (That is, for those who already know, Rq is 
one half the Schwartzschild radius of our sun.) In this equation p and e still carry dimensions 
of, say, ergs/cm^. Therefore, let us define dimensionless energy density and pressure, e and 
P, by 

p = eop, e = eoe (19) 

where eo has dimensions of energy density. This eo is not the same as defined in Eq. (11). 
Its numerical choice here is arbitrary, and a suitable strategy is to make that choice based 
on the dimensionful numbers that define the problem at hand. We'll employ this strategy 
to fix it below. For a polytrope, we can write 

p = Ke'^ , where K = Ke^^^ is dimensionless. (20) 

It is easier to solve Eq. (|T8|l for p, so we should express e in terms of it, 

-e={p/Kf'^. (21) 

Equation p8|) can now be recast in the form 

dp{r) ^ ap{rY/^'M{r) 
dr 

where the constant a is 

a = R^/K'/'' = Ro/{KerY"'. (23) 

Equation ()22|1 has dimensions of 1/km, with a in km (since Rq is). That is, it is to be 
integrated with respect to r, with r also in km. 

We can choose any convenient value for a since eo is still free. For a given value of a, eo 
is then fixed at 

-,1-7 

€0^ -{ — ] ■ (24) 
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We also need to cast the other coupled equation, Eq. in terms of dimensionless p and 
M, 

'-^-PrMrY\ (25) 

where 

^ 47reo ^ Airep 
^ M0c2i?V7 M0c2(i^e2"')V7 ■ ^ 

Equation (j^^ also carries dimensions of 1/km, the constant (3 having dimesnions 1/km^. 
Note that, in integrating out from r = 0, the initial value of A4{0) = 0. 



D. Integrating the Polytrope Numerically 

Our task is to integrate the coupled first-order differential equations (DE), Eqs. (j22|l and 
I^Kll . out from the origin, r = 0, to the point R where the pressure falls to zero, p{R) = 
18l |. To do this we need two initial values, p{0) (which must be positive) and Ai{0) ( which 
we already know must be 0). The star's radius, R, and its mass M = Ai{R) in units of Mq 
will vary, depending on the choice for p{0). 

For purposes of numerical stability in solving Eqs. ()22p and ()25p . we want the constants 
a and (3 to be not much different from each other (and probably not much different from 1). 
We will see below that this can be arranged for both of the two polytropic EoS's discussed 
above for white dwarfs. 

Our coupled DE's are quite non-linear. In fact, because of the p^^"' factors, the solution 
will become complex when p{r) < 0, i.e., when r > R. Thus we will want to recognize when 
this happens. How can this be programmed? 

Mathematica and similar symbolic/numerical packages have built-in first-order DE 
solvers. Perhaps the solver is as simple as a fixed, equal-step Runge-Kutta routine (as 
in MathCad 7 Standard), but there are often more sophisticated solvers in more recent ver- 
sions. These packages also allow for program control constructs such as do-loops, whiles and 
the like. 

Thus, consider a do-loop on a variable f running in appropriately small steps over a 
range that is sure to contain the expected value of R. Call the DE solver inside this loop, 
integrating the coupled DE's from r = to f. When the solver routine exits, check to see if 
the last value of p, i.e., p(f), has a real part which has gone negative. If so, then break out 
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of the loop, calling R = f. If not, go on to the next larger value of f and call the DE solver 
again. 

More discussion of how to program the integration of the DE's is inappropriate here, 
since we want to encourage the student to learn from programming to appreciate how the 
symbolic/numerical package is used. 

E. The Relativistic Case kp ^ rUg 

This is the regime for white dwarfs with the largest mass. A larger mass needs a greater 
central pressure to support it. However, large central pressures mean the squeezed electrons 
become relativistic. 

Recall that the polytrope exponent 7 = 4/3 for this case and the equation of state is 
given by P = Kj.f.ie'^ with Kj-^i given by Eq. After some trial and error, we chose in our 
program (the student may want to try something else) 

a = Ro = 1.473 km [kp > m^] , (27) 

which in turn fixes, from Eq. (|23), 

eo = 7.463 x lO^^ergs/cm^ = 4.17 MgcVkm^ [kp > me] . (28) 

The first question the student should ask, in checking this number, is whether such a large 
number is physically reasonable. 

Continuing with the kp ^ rrie numerics, Eqs. (fT^ and give 

P = 52.46 /km^ [kp > m^] , (29) 

which is about 30 times larger than a, but probably manageable from the standpoint of 
performing the numerical integration. 

In our first attempt to integrate the coupled DE's for this case (using a do-loop as 
described above) we chose p{0) = 1.0. This gives us a white dwarf of radius R ^ 2 km, 
which is miniscule compared with the expected radius of ~ 10'^ km! Why? What went 
wrong? 

The student who also makes this kind of mistake will eventually realize that our choice 
of scale, eo = 4.17M0C^/km^, represents a huge energy density. One can simply estimate the 
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TABLE I: Radius R (in km) and mass M (in Mq) for white dwarfs with a relativistic electron 
Fermi gas EoS. 



p{0) 


R 


M 


10-14 


4840 


1.2431 


10-15 


8600 


1.2432 


10-16 


15080 


1.2430 



average energy density of a star with a 10'^ km radius and a mass of one solar mass by the 
ratio of its rest mass energy to its volume, 

(e) ^ ^ = 10-^^ MocVW , (30) 

which is much, much smaller than the eo here. In addition, the pressure p is about 2000 
times smaller than the energy density e (see Fig. Thus, choosing a starting value of 
p{0) ~ 10~i5 would probably be more physical. 

Doing so does give much more reasonable results. Table H] shows our program's results 
for R and M and how they depend on p{0). The surprise here is that, within the numerical 
error expected, all these cases have the same mass! Increasing the central pressure doesn't 
allow the star to be more massive, just more compact. 

It turns out that this result is correct, i.e., that the white dwarf mass is independent of 
the choice of the central pressure. It is not easy to see this, however, from the numerical 
integration we have done here. The discussion in terms of Lane-Emden functions Q| shows 
why, though the mathematics here might be a bit steep for many undergraduates. For this 
reason, we quote without proof the analytic results. For the case of a polytropic equation 
of state p = Ke^ , the mass 

and the radius 

In the above-mentioned solutions, Ci and ^(Ci) ^'^^ numerical constants that depend on the 
polytropic index 7. By examining Eq. (j3ip . we see that for 7 = 4/3 the mass is independent 
of the central energy density, and hence also the central pressure Pq. Also, note that from 
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FIG. 3: p{r) and 7W(r) for white dwarfs using the relativistic electron Fermi gas model. Here 
p{0) = 10-16. 

Eq. (jH^ . the radius decreases with increasing central pressure as oc p^^ = Pq^^'^- 

In any case, the student should notice this point and use it as check of the numerical 
results obtained. Figure El shows the dependence of p{r) and A4{r) on distance for the case 
p{0) = IQ-i^. It is interesting that p{r) becomes small and essentially fiat around 8000 km 
before finally going through zero at i? = 15, 080 km. 

The results and graphs shown here were generated with a Mathematica 4.0 program, but 
we were able to reproduce them using MathCad 7 Standard. In that case, however, program- 
ming a loop is difficult, so we searched by hand for the endpoint (where the real part of p{r) 
goes negative). More recent versions of MathCad have more complete program constructs, 
such as while-loops, so this process could undoubtedly be automated. (Alternatively, the 
student might try to solve for a root of p{r) = 0.) 
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F. The Non-Relativistic Case, kp <^ nie 



Eventually, as the central pressure p{0) gets smaller, the electron gas is no longer rela- 
tivistic. Also as the pressure gets smaller, it can support less mass. This moves us in the 
direction of the less massive white dwarfs, and, as it turns out, these dwarfs are larger (in 
radius) than the ones in the last section. 

In the extreme case, when kp <^ mg, we can integrate the structure equations for the 
other polytropic EoS, where 7 = 5/3. The programming for this is very much the same as 
in the 4/3 case, but the numbers involved are quite different (as are the results). 

Inserting the values of the physical constants in Eq. (|17|). we find 

90 cm^ , , 

ir„on-rel = 3.309 X lO'^^^ ^ . 33 

ergs /'^ 

This time, however, and after some experimentation, we chose the constant 

a = 0.05 km [kp < nie] , (34) 

which then fixes 

eo = 2.488 x lO^Vgs/cm^ = 0.01392 MqC^ /km^ [kp < m^] . (35) 

Note that this eo is much smaller than our choice for the relativistic case. The other constant 
we need, from Eq. (j^ . is 

p = 0.005924 /km^ [kp < rrie] , (36) 

which, unlike the relativistic case, is not larger than a but smaller. 

When we first ran our Mathematica code for this case, we (inadvertantly) tried a value 
of p{0) = 10^^^. This gave us a star with radius R = 5310 km and mass M = 3.131. Oops!, 
that mass is bigger than the largest mass of 1.243 that we found for the relativistic EoS! 
What did we do wrong? 

What happened (and the student can set up her program so this trap can be avoided) 
is that the choice p{0) = 10~^^ violates the assumption that kp <^ nif,. One really needs 
values for p(0) < 4 x 10^^^. This says, in fact, that the relativistic p(0) = 10~^^ case that 
we plotted in Fig. |21is not really relativistic. 

Results for the non-relativistic case for the last two values of p(0) in Table H] are shown in 
Table ITU It is quite instructive to compare the differences in the two tables. The masses are, 

16 



TABLE II: Radius R (in km) and mass M (in Mq) for white dwarfs with a non-relativistic electron 
Fermi gas EoS. 



p{0) 


R 


M 


10-15 


10620 


0.3941 


10-16 


13360 


0.1974 




5000 10000 



FIG. 4: p{r) for a white dwarf using the non-relativistic electron Fermi gas model with central 
pressure p{0) = lO^^^. 

of course, smaller, as expected, and now they vary with p{0). Somewhat surprising is that 
the non-relativistic radius is bigger for p{0) = IQ-^^ smaller for p{0) = 10~^^. Figured 
shows the pressure distribution for the latter case, to be compared with the corresponding 
graph in Fig. IHl Note that this pressure curve does not have the peculiar, long flat tail 
found using the relativistic EoS. 

In fact, by this time the student should have realized that neither of these polytropes is 
very physical, at least not for all cases. The non-relativistic assumption certainly does not 
work for central pressures p{0) > 10~^^, i.e., for the more massive (and more common) white 
dwarfs. On the other hand, the relativistic EoS certainly should not work when the pressure 
becomes small, i.e., in the outer regions of the star (where it eventually goes to zero at the 
star's radius). So, can one find an EoS to cover the whole range of pressures? 

We haven't actually done this for white dwarfs, but the program would be much like that 
discussed below for the full neutron star. Given the transcendental expressions for energy 
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and pressure that generate the curve shown in Fig. it should be possible to find a fit 
(using, for example, the built-in fitting function of Mathematica) like 

-e{p) = A^^p'/' + A^p'/' . (37) 

The second term dominates at high pressures (the relativistic case), but the first term takes 
over for low pressures when the kp ^ trif. assumption does not hold. (Setting the two terms 
equal and solving for p, as Chandrasekhar and Fowler did, gives the value of p when special 
relativity starts to be important.) This expression for e{p) could then be used in place of the 
p^/"^ factors on the right hand sides of the structure equations. Proceed to solve numerically 
as before. We leave this as an exercise for the interested student. 

IV. PURE NEUTRON STAR, FERMI GAS EOS 

Having by now become warm, the student can now tackle neutron stars. Here one must 
include the general relativistic (GR) contributions represented by the three dimensionless 
factors in the TOV equation, Eq. One of the first things that comes to mind is how one 
deals numerically with the (apparent) singularities in these factors at r = 

Also, as in the case of the white dwarfs, there is a question of what to use for the EoS. 
In this section we show what can be done for pure neutron stars, once again using a Fermi 
gas model for, now, a neutron gas instead of an electron gas. Such a model, however, is 
unrealistic for two reasons. First, a real neutron star consists not just of neutrons but 
contains a small fraction of protons and electrons (to inhibit the neutrons from decaying 
into protons and electrons by their weak interactions). Second, the Fermi gas model ignores 
the strong nucleon-nucleon interactions, which give important contributions to the energy 
density. Each of these points will be dealt with in sections below. 

A. The Non-Relativistic Case, kp <^ mn 

For a pure neutron star Fermi gas EoS one can proceed much as in the white dwarf case, 
substituting the neutron mass m„ for the electron mass rrie in the equations found in Sec. 
3. When kp <^ m„ one finds, again, a polytrope with 7 = 5/3. (More exercises for the 
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TABLE III: Radius R (in km) and mass M (in Mq) for pure neutron stars with a non-relativistic 
Fermi gas EoS. 



p{0) 


R (Newton) 


M (Newton) 


R (GR) 


M (GR) 


10-4 


16.5 


0.7747 


15.25 


0.6026 


10"^ 


20.8 


0.3881 


20.00 


0.3495 


lo-'' 


26.3 


0.1944 


25.75 


0.1864 



student.) The K corresponding to that in Eq. (fTTjl is 

^2 / o„2r7 \ 5/3 



3n'Z 



non— rel 



6.483 X 10 



-26 



cm 



ergs 



2/3 ■ 



This time, choosing a = 1 km, one finds the scahng factor from Eq. ()24|1 to be 

eo = 1.603 X 10^^ ergs/cm^ = 0.08969 MocVkm^ ■ 
Further, from Eqs. (j2I| and fl^ . 

K = 1.914 and P = 0.7636 /km^ . 



(38) 



(39) 



(40) 



Note that, in this case, the constants a and (3 are of similar size. 

Making an estimate of the average energy density of a typical neutron star (mass = Mq, 
R = 10 km), one expects that a good starting value for the central pressure p{0) to be of 
order 10-^ or less. Our program for this situation is essentially the same as the one for 
non-relativistic white dwarfs but with appropriate changes of the distance scale. It gives the 
results shown in Table lllll Note that the GR effects are small, but not negligible, for this 
non-relativistic EoS. As in the white dwarf case, these are the smaller mass stars. One sees 
that as the mass gets smaller, the gravitational attraction is less and thus the star extends 
out to larger radii. 



B. The Relativistic Case, kp » m„ 

Here there is again a polytropic EoS, but with 7 = 1. In fact, p = e/3, a well-known 
result for a relativistic gas. The conversion to dimensionless quantities becomes very simple 
in this case with relationships like K = K = 1/3. It is still useful to factor out an eo, which 
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in our program we took to have a value 1.6 x 10^^ erg/cm^, as suggested by the value in the 
previous sub-section. Then, if we choose this time 



a = 3Ri 



-0 — 



4.428 km 



(41) 



we find 



P = 3.374 /km^ . 



(42) 



We expect central pressures p{0) in this case to be greater than 10 ^. Other than these 
changes, we wrote a similar program to the one above, taking care to avoid exponents like 



Running that code gives, at first glance, enormous radii, values of R greater than 50 
km! We can imagine the student looking frantically for a program bug that isn't there. In 
fact, what really happens is that, for this EoS, the loop on f runs through its whole range, 
since the pressure p{r) never passes through zero. (A plot of p{r) looks quite similar, but 
for distance scale, to that shown in Fig. IHl where 7 = 4/5.) It only falls monotonically 
toward zero, getting ever smaller. By the time the student recognizes this, she will probably 
also have realized that the relativistic gas EoS is inappropriate for such small pressures. 
Something better should be done (as in the next sub-section). 

It turns out that the structure equations for 7 = 1 are sufficiently simple that an analytic 
solution for p{r) can be found, which corroborates the above remarks about not having a 
zero at a finite R. A suggestion for the student is to try a power-law Ansatz. 

C. The Fermi Gas EoS for Arbitrary Relativity 

In order to avoid the trap of the relativistic gas, one should find an EoS for the non- 
interacting neutron Fermi gas which works for all values of the relativity parameter x = 
kp/iTLnC. Taking a hint from the two polytropes, one can try to fit the energy density as a 
function of pressure, each given as a transcendental function of kp, with the form 



For low pressures the non-relativistic first term dominates over the second. (The power in 
the relativistic term is changed from that in Eq. (p?7jl .) It is again useful to factor out an eo 



1/(7-1). 



e(p) = At^Rp^'^ + Arp . 



(43) 
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from both e and p. In this case, it is more natural to define it as 

mtc^ n(. ergs Mp^c^ , . 

eo = " = 5.346 x 10^^ —2- = 0.003006 . (44 

yi-K^Tiy^ cm'^ km 

Mathematica can easily create a table of exact e and p values as a function oikp. The 
dimensionless A-values can then be fit using its built-in fitting function. From our efforts 
we found, to an accuracy of better than 1% over most of the range of hp |2o| . 

Anr = 2.4216 , = 2.8663 . (45) 

We used the fitted functional form for e of Eq. (j43p in a Mathematica program similar 
to that for the neutron star based on the non-relativistic EoS. With the eo of Eq. ()44|1 and 
choosing a = = 1.476 km, we obtain /3 = 0.03778. Our results for a starting value of 
p(0) = 0.01, clearly in the relativistic regime, are 

R = 15.0, M = 1.037, Newtonian equations (46) 
i? = 13.4 , M = 0.717, full TOV equation . (47) 

For this more massive star, the GR effects are significant (as should be expected from the 
size of GM/c^R, about 10% in this case). Figure El displays the differences. 

It is now instructive to make a long run of calculations for a range of p{0) values. We 
display in Fig.lHla (parametric) plot of M and R as they depend on the central pressure. The 
low-mass/large-radius stars are to the right in the graph and correspond to small starting 
values of p{0). As the central pressure increases, the total mass that the star can support 
increases. And, the bigger the star mass, the bigger the gravitational attraction which 
draws in the periphery of the star, making stars with smaller radii. That is, increasing p{0) 
corresponds to "climbing the hill," moving upward and to the left in the diagram. 

At about p{0) = 0.03 one gets to the top of the hill, achieving a maximum mass of about 
0.8 Mq at a radius R ^ 11 km. That maximum M and its R agree with Oppenheimer and 
Volkov's seminal 1939 result for a Fermi gas EoS. 

What about the solutions in Fig. IHlthat are "over the hill," i.e., to the left of the maxi- 
mum? It turns out that these stars are unstable against gravitational collapse into a black 
hole. The question of stability, however, is a complicated issue [2^, perhaps too difficult for 
a student at this level. The fact that things collapse to the left of the maximum, however, 
means that one probably shouldn't worry too much about the peculiar curly-cue tail to the 
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FIG. 5: p{r) and A^(r) (r in km) for a pure neutron star with central pressure p(0) = 0.01, using a 
Fermi gas EoS fit valid for all values ofkp- The thin curves are results from the classical Newtonian 
structure equations, while the thick ones include general relativistic corrections. 

M-R curve in the figure. It appears to be an artifact for very large values of p(0), also seen 
in other calculations, even though it is especially prominent for this Fermi gas EoS. 

D. Why Is There a Maximum Mass? 

On general grounds one can argue that cold compact objects such as white dwarfs and 
neutron stars must possess a limiting mass beyond which stable hydrostatic configurations 
are not possible. This limiting mass is often called the maximum mass of the object and 
was briefly mentioned in the discussion at the end of Sec. 2.2 and that relating to Fig. 6. In 
what follows, we outline the general argument. 
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M{R) 

Ir 

0.8 




0.2 

5 10 15 20 25 

FIG. 6: The mass M (in Mq) and radius R (in km) for pure neutron stars, using a Fermi gas 
EoS. The stars of low mass and large radius are solutions of the TOV equations for small values 
of central pressure p{0)- The stars to the right of the maximum aX R = 11 are stable, while those 
to the left will suffer gravitational collapse. 

The thermal component of the pressure in cold stars is by definition negligible. Thus, 
variations in both the energy density and pressure are only caused by changes in the density. 
Given this simple observation, let us examine why we expect a maximum mass in the 
Newtonian case. 

Here, an increase in the density results in a proportional increase in the energy density. 
This results in a corresponding increase in the gravitational attraction. To balance this, 
we require that the increment in pressure is large enough. However, the rate of change 
of pressure with respect to energy density is related to the speed of sound (see Sec. 6.3). 
In a purely Newtonian world, this is in principle unbounded. However, the speed of all 
propagating signals cannot exceed the speed of light. This then puts a bound on the pressure 
increment associated with changes in density. 

Once we accept this bound, we can safely conclude that all cold compact objects will 
eventually run into the situation in which any increase in density will result in an additional 
gravitational attraction that cannot be compensated for by the corresponding increment in 
pressure. This leads naturally to the existence of a limiting mass for the star. 

When we include general relativistic corrections, as discussed in Sec. 2.2 earlier, they act 
to "amplify" gravity. Thus we can expect the maximum mass to occur at a somewhat lower 
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mass than in the Newtonian case. 



V. NEUTRON STARS WITH PROTONS AND ELECTRONS, FERMI GAS EOS 

As mentioned at the beginning of the last section, neutron stars are not made only of 
neutrons. There must also be a small fraction of protons and electrons present. The reason 
for this is that a free neutron will undergo a weak decay, 

n ^ p + e~ + Ue , (48) 

with a lifetime of about 15 minutes. So, there must be something that prevents this decay 
in the case of the star, and that is the presence of the protons and electrons. 

The decay products have low energies (m„ — rrip — mg = 0.778 MeV), with most of that 
energy being carried away by the light electron and (nearly massless) neutrino [2^. If all 
the available low-energy levels for the decay proton are already filled by the protons already 
present, then the Pauli exclusion principle takes over and prevents the decay from taking 
place. 

The same might be said about the presence of the electrons, but in any case the electrons 
must be present within the star to cancel the positive charge of the protons. A neutron star 
is electrically neutral. We saw earlier that the number density of a particle species is fixed 
in terms of that particle's Fermi momentum [see Eq. ((7j)]. Thus equal numbers of electrons 
and protons implies that 

kF,p = kF,e ■ (49) 

In addition to charge neutrality, we also require weak interaction equilibrium, i.e., as 
many neutron decays [Eq. (jlHjl ] taking place as electron capture reactions, p + e~ ^ n + Ug. 
This equilibrium can be expressed in terms of the chemical potentials for the three particle 
species, 

IJ'n = fJ'p + f^e ■ (50) 

We already defined the chemical equilibrium for a particle in Sec. 3.2 after Eq. (fT^ . 

l^i{kF,i) = ^ = i^h + ^ly^^ , i = n,p,e. (51) 

where, for the time being, we have set c = 1 to simplify the equations somewhat. (The 
student is urged to prove the right-hand equality.) From Eqs. (|^. and we can 
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find a constraint determining kp^p for a given kp^n, 

(4n + ^ly^' - (4. + ^ly^" - (4. + miY^' = . (52) 

While an ambitious algebraist can probably solve this equation for kp^p as a function of kp^n, 
we were somewhat lazy and let Mathematica do it, finding 

kFAkF,n) = 2(fc^,„ + OV2 (53) 

kp„ +m^-m^ „ m. 



The total energy density is the sum of the individual energy densities, 

etot = X! ' (^^) 



i=n,p,e 



where 



eiikp,i) = ' (fc^ + miY''^k'^dk = eo ei(xi, , (56) 
Jo 

eo = m^/STT^/i^ (57) 

Mx,,y,) = r\u^ + y^if'Wdu, (58) 

JO 

a^i = kp^i/rrii, yi = mi/mn. (59) 
The corresponding total pressure is 

Ptot = ^ Pi, (60) 



and, as before 



0, 



i=n,p,e 
rkf 



Pi 



Pi{kp,i) = / {k"^ + mi)~^^'^k^dk = eoPiixi,yi) , (61) 
Jo 

i{xi,yi) = r\u^ + y^)-'/\^du. (62) 
Jo 

Using Mathematica the (dimensionless) integrals can be expressed in terms of log and 
sinh~^ functions of Xj and yi. One can then generate a table of etot versus ptot values for an 
appropriate range of kp^nS. This, in turn, can be fitted to the same form of two terms as 
before in Eq. pHjl . We found, this time, the coefficients to be 

Anr = 2.572 , Ak = 2.891 . (63) 

These coefficients are not much changed from those in Eq. for the pure neutron star. 
Therefore, we expect that the M versus R diagram for this more realistic Fermi gas model 
would not be much different from that in Fig. IHl 
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VI. INTRODUCING NUCLEAR INTERACTIONS 



Nucleon-nucleon interactions can be included in the EoS (they are important) by con- 
structing a simple model for the nuclear potential that reproduces the general features of 
(normal) nuclear matter. In doing so we were much guided by the lectures of Prakash . 

We will use MeV and fm (10^^^ cm) as our energy and distance units for much of this 
section, converting back to Mq and km later. We will also continue setting c = 1 for now. 
In this regard, the important number to remember for making conversions is he = 197.3 
MeV-fm. We will also neglect the mass difference between protons and neutrons, labeling 
their masses as m^r. 

The von Weizacker mass formula 2J| for nuclides with Z protons and neutrons gives. 



for normal symmetric nuclear matter {A = N + Z with N = Z), an equilibrium number 
density uq of 0.16 nucleons/fm^. For this value of no the Fermi momentum is kp = 263 
MeV/c [see Eq. ((7j)]. This momentum is small enough compared with rriN- = 939 MeV/c^ 
to allow a non-relativistic treatment of normal nuclear matter. At this density, the average 
binding energy per nucleon, BE = E / A—itln^ is —16 MeV. These are two physical quantities 
we definitely want our nuclear potential to respect, but there are two more that we'll need 
to fix the parameters of the model. 

We chose one of these as the nuclear compressibility, Kq, to be defined below. This is a 
quantity which is not all that well established but is in the range of 200 to 400 MeV. The 
other is the so-called symmetry energy term, which, when Z = 0, contributes about 30 MeV 
of energy above the symmetric matter minimum at uq. (This quantity might really be better 
described as an asymmetry parameter, since it accounts for the energy that comes in when 

Ny^Z.) 



A. Symmetric Nuclear Matter 

We defer the case when N ^ Z, which is our main interest in this paper, to the next 
sub-section. Here we concentrate on getting a good (enough) EoS for nuclear matter when 
N = Z, or, equivalently, when the proton and neutron number densities are equal, n„ = Up. 
The total nucleon density n = nn + Up. 

We need to relate the first three nuclear quantities, no, BE, and Kq to the energy density 
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for symmetric nuclear matter, e(n). Here n = nlkp) is the nuclear density (at and away 
from uq). We will not worry in this section about the electrons that are present, since, as 
was seen in the last section, its contribution is small. The energy density now will include 
the nuclear potential, V{n), which we will model below in terms of two simple functions 
with three parameters that are fitted to reproduce the above three nuclear quantities. [The 
fourth quantity, the symmetry energy, will be used in the next sub-section to fix a term in 
the potential which is proportional to (A^ — Z)/A.] 

First, the average energy per nucleon, E/A, for symmetric nuclear matter is related to e 

by 

E{n)/A = e{n)/n, (64) 

which includes the rest mass energy, rriN and has dimensions of MeV. As a function of n, 
E{n)/A — m^v has a minimum at n = no with a depth BE = —16 MeV. This minimum 
occurs when 

d f E(n)\ d f e{n)\ ^ 

' ^ ^ > - ' ^ ^ ' - at n = no . (65) 



dn \ A J dn \ n ^ 
This is one constraint of the parameters of V{n). Another, of course, is the binding energy, 

e(n) 

mj^ = BE at n = hq . (66) 

n 

The positive^ curvature at the bottom of this valley is related to the nuclear (in)compress- 
ibility by 



K{n) = 9^ = 9 
dn 



{-] + 2n 



d' /e\ d 



dii? \nj ' dn \nj ^ ^"^ 
using Eq. ()13|1 . which defines the pressure in terms of the energy density. At n = no this 
quantity equals Kq. (The factor of 9 is a historical artifact from the convention originally 
defining Kq.) 

(Question: why does one not have to calculate the pressure at n, = no?) 
The N = Z potential in e(n) we will model as 

e(n) 3h'^kl A B ^ 

n 5 2mjv 2 cr + 1 

where u = n/n^ and a are dimensionless and A and B have dimensions of MeV. The first 
term represents the rest mass energy and the second the average kinetic energy per nucleon. 
[These two terms are leading in the non-relativistic limit of the nucleonic version of Eq. 
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()1U|).] For kpino) = k% we will abbreviate the kinetic energy term as (-E"^), which evaluates 
to 22.1 MeV. The kinetic energy term in Eq. (jUHj) can be better written as {E'p) v?^'^. 

From the above three constraints, Eqs. (jU^ - (jH7|) . and noting that m = 1 at n = no, we 
get three equations for the parameters A, B, and a: 

B 



E%) + --^ + 



z,{E%) + '-^ + 



A 

2" 

A 

2" 



a + 1 
Ba 



BE, 



0, 



9 



a + 1 

A + Ba = ^ 
9 



(69) 
(70) 
(71) 



Solving these equations (which we found easier to do by hand than with Mathematica) , one 
finds 



(72) 

(73) 
(74) 



a = 


Ko + 2 ) 




-9BE ' 


B = 


a + l 
a-1 




A = 


BE - 





-BE 
B. 

Numerically, for Kq = 400 MeV (which is perhaps a high value), 

A = -122.2 MeV, B = 65.39 MeV, a = 2.112 . 



(75) 



Note that a > 1, a point we will come back to below, since it violates a basic principle of 
physics called "causality." 

The student can try other values of Kq to see how the parameters A, B, and a change. 
More interesting is to see how the interplay between the A- and B-terms gives the valley at 
n = riQ. Figure [7| shows E / A — ttln as a function of n using the parameters of Eq. (f73j) . We 
would hope the student notices the funny little positive bump in this plot near n = and 
sorts out the reason for its occurrence. 

Given e(n) from Eq. one can find the pressure using Eq. (fT!^ . 



pin) = n^— (— 
an \n 



no 



u 



5/3 



A 



B 



a 



+ 77« + 

2 a + l 



.cr+l 



U 



(76) 



For the parameters of Eq. (ffH|l its dependence on n is shown in Fig. |H1 On first seeing this 
graph, the student should wonder why p{u = 1) = p(no) = 0. And, what is the meaning of 
the negative values for pressure below u = 11 (Hint: what is "cavitation"?) 
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EjA - ruN 




FIG. 7: The average energy per nucleon, less its rest mass, as a function of n = n/ng (in MeV). 
The position of the minimum is at n = no = 0.16 fm~^, its depth there is BE = — 16 MeV, and 
its curvature (second derivative) there corresponds to nuclear compressibility Kq = 400 MeV. 

p{u) 




FIG. 8: The pressure for symmetric nuclear matter as a function of u = n/riQ. The student should 
ask what it means when the pressure is negative and why it is at n = 1. 

So, if this N = Z case were all we had for the nuclear EoS, a plot of e{n) versus p{n) 
would only make sense for n > uq. Such a plot looks much like a parabola opening to the 
right for the range < m < 3. At very large values of m, however, e ~ p/3, as it should for 
a relativistic nucleon gas (see Section 4.2). We don't pursue this symmetric nuclear matter 
EoS further since our interest is in the case when N ^ Z 
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B. Non-Symmetric Nuclear Matter 

We continue following Prakash's notes closely. Let us represent the neutron and proton 
densities in terms of a parameter a as 



1 + a 
2 ^ 



1 — a 



■ n . 



This a is not to be confused with the constant defined in Eq. 
a = 1. Note that 



(77) 

For pure neutron matter 

A ' 

so we can expect that the isospin-symmetry-breaking interaction is proportional to a (or 
some power of it). An alternative notation is in terms of the fraction of protons in the star, 

1 — a 



a 



Un — Up N — Z 



rir. 



X 



(79) 



n 2 

We now consider how the energy density changes from the symmetric case discussed above, 
where a = (or x = 1/2). 

First, there are contributions to the kinetic energy part of e from both neutrons and 
protons. 



eKE[n,a) 



3 ^F,n 3 k"^ 



5 2m 



TV 



I F,p 

5 zniN 



:i + a 



,5/3 



'1 — a 



,5/3 



where 



(Ef) = - 



3 fSn'nV^' 



(80) 



(81) 



5 2mN \ 2 J 

is the mean kinetic energy of symmetric nuclear matter at density n. For n = no we note 
that {Ep) = 3 {Ep) /5 [see Eq. For non-symmetric matter, a 7^ 0, the excess kinetic 

energy is 



AeKE{n,a) = eKsin, a) - eKE{n,0) 



[l + a 



,5/3 



'1 — a 



,5/3 



n ' 



(Ef) {2 

For pure neutron matter, a = 1, 

AeKEin,a 



2/3 



(l_a;)5/3^^5/3] 



n 



(82) 



(83) 
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It is also useful to expand to leading order in a, 

AeKE{n,a) = n{Ep)-a^ il + — + ■■■] (84) 

Keeping terms to order is evidently good enough for most purposes. For pure neutron 
matter, the energy per particle (which, recall, is e/n) at normal density is AexEino, l)/'no ~ 
13 MeV, more than a third of the total bulk symmetry energy of 30 MeV, our fourth nuclear 
parameter. 

Thus the potential energy contribution to the bulk symmetry energy must be 20 MeV or 
so. Let us assume the quadratic approximation in a also works well enough for this potential 
contribution and write the total energy per particle as 

E{n,a) = E{n,0) + a^S{n), (86) 

The isospin-symmetry breaking is proportional to a^, which reflects (roughly) the pair- wise 
nature of the nuclear interactions. 

We will assume S{u), u = n/uQ, has the form 

Siu) = (2^/3 - 1)^ ) {u'/' - Fiu)) + SoFiu) . (87) 

Here 5*0 = 30 MeV is the bulk symmetry energy parameter. The function F{u) must satisfy 
F{1) = 1 [so that S{u = 1) = 5*0] and -F(O) = [so that S{u = 0) = 0; no matter means 
no energy]. Besides these two constraints there is, from what we presently know, a lot of 
freedom in what one chooses for F{u). We will make the simplest possible choice here, 
namely, 

F{u) = u , (88) 

but we encourage the student to try other forms satisfying the conditions on F{u), such as 
a/u, to see what difference it makes. 

Figure IHl shows the energy per particle for pure neutron matter, E{n,l) — rriN, as a 
function of u for the parameters of Eq. ()75|1 and 5*0 = 30 MeV. In contrast with the a = 
plot in Fig. E{n,l) > and is monotonically increasing. The plot looks almost quadratic 
as a function of u. The dominant term at large u goes like m"", and a = 2.112 (for this case). 
However, one might have expected a linear increase instead. We will return to this point in 
Sec. 6.3. 
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FIG. 9: The average energy per neutron (less its rest mass), in MeV, for pure neutron matter, as 
a function of n = n/riQ. The parameters for this curve are for a nuclear compressibility Kq of 400 
MeV. 

Given the energy density, e(n, a) = nouE{n,a), the corresponding pressure is, from Eq. 

da, 



p(n,x) = u—e(n,a) — e(n,a) 
du 



p{n, 0) + noa^ 



22/3 



- ) {2u'/' - 3u') + Sou' 



(89) 



where p{n,0) is defined by Eq. (f7H|) . Figure ITUl shows the dependence of the pure neutron 
p{n, 1) and e(n, 1) on u = n/rio, ranging from to 10 times normal nuclear density. Both 
functions increase smoothly and monotonically from u = 0. We hope the student would 
wonder why the pressure becomes greater than the energy density around u = 6. Why 
doesn't it go like a relativistic nucleon gas, p = e/3? (Hint: check the assumptions.) 

One can now look at the EoS, i.e., the dependance of p on e (the points in Fig. [TT|). The 
pressure is smooth, non-negative, and monotonically increasing as a function of e. In fact it 
looks almost quadratic over this energy range (0 < m < 5). This suggests that it might be 
reasonable to see if one can make a simple, polytropic fit. If we try that using a form 



we find the fit shown in Fig. ^2 as the solid curve with 

= 3.548 X 10-^ 7 = 2.1 
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FIG. 10: The pressure (dashed curve) and energy density (sohd) for pure neutron matter, as a 
function of u = n/riQ. Units for the y-axis are MeV/fm^. This curve uses parameters based on a 
nuclear compressibihty Kq = 400 MeV. 

where kq has appropriate units so that p and e are in MeV/fm^. (We simply guessed and 
set 7 to that value.) 

This polytrope can now be used in solving the TOV equation for a pure neutron star with 
nuclear interactions. Alternatively, one might solve for the structure by using the functional 
forms from Eq. (jHEI), multiplied by n, and Eq. ()89j) directly. We defer that for a bit, since 
it would be a good idea to first find an EoS which doesn't violate causality, a basic tenet of 
special relativity. 



C. Does the Speed of Sound Exceed That of Light? 

What is the speed of sound in nuclear matter? Starting from the elementary formula for 
the square of the speed of sound in terms of the bulk modulus j27|, one can show that 

/c,\2 ^B_^dp^ dp/dn . 
\ c J pc^ de de / dn 

To satisfy relativistic causality we must require that the sound speed does not exceed that 
of light. This can happen when the density becomes very large, i.e., when u oo. For the 
simple model of nuclear interactions presented in the last section, the dominant terms at 
large u in p and e are those going like u'^^'^ . Thus, from Eq. ()8(i|l . multiplied by n, and Eq. 
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FIG. 11: The equation of state for pure neutron matter (a = 1), i.e., the dependence of pressure 
(y-axis) versus energy density (x-axis). Units for both axes are MeV/fm^, and the nuclear com- 
pressibiUty in this case is Kq = 400 MeV. The points are values calculated directly from Eq. ()86() . 
multiplied by n, and Eq. (|89|) . while the solid curve is a fit to these points given in Eqs. (|90|) and 
(|FT|1. 

we see that 

'^^ = ^^^ = 2.11 (93) 

c / dejdn 

for the parameters of Eq. (jTS)), and indeed for any set of parameters with Kq greater than 
about 180 MeV. 

One can recover causality (i.e., speeds of sound less than light) by assuring that both 
t{u) and p{u) grow no faster than v?. There must still be an interplay between the A- and 
S-terms in the nuclear potential, but one simple way of doing this is to modify the 5-term 
by introducing a fourth parameter C so that, for symmetric nuclear matter (a = 0), 

W",0) = ^. + ^^^f^. (94) 

One can choose C small enough so that the effect of the denominator only becomes ap- 
preciable for very large u. The presence of the denominator would modify and complicate 
the constraint equations for A, S, and a from those given in Eqs. ()(i9|l - (f7T|l . However, for 
small C, which can be chosen as one wishes, the values for the other parameters should not 
be much changed from those, say, in Eq. (j75|) . Thus, with a little bit of trial and error, 
one can simply readjust the A, B, and a values to put the minimum of E/ A — rriN dX the 
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right position {uq) and depth (BE), hoping that the resulting value of the (poorly known) 
compressibility Ko remains sensible. 

In our calculations we chose C = 0.2 and started the hand search with the Kq = 400 
MeV parameters in Eq. (j75|) . We found that, by fiddling only with B and a, we could re-fit 
hq and B with only small changes, 

B = 65.39 ^ 83.8 MeV , a = 2.11 ^ 2.37 , (95) 

somewhat larger than before. For these new values of B and a, A changes from -122.2 MeV 
to -136.7 MeV, and Kq from 400 to 363.2 MeV. That is, it remains a reasonable nuclear 
model. 

One can now proceed as in the last section to get e(n, a), p{n,a), and the EoS, p{e,a). 
The results are not much different from those shown in the figures of the previous sub- 
section. This time we decided to live with a quadratic fit for the EoS for pure neutron 
matter, finding 

p(e, 1) = Koe^ , Ko = 4.012 x 10"^ . (96) 

This is not much different from before, Eq. 1)9111 . Somewhat more useful for solving the TOV 
equation is to express e in terms of p, 

e{p) = {p/Ko)f' . (97) 



D. Pure Neutron Star with Nuclear Interactions 



Having laid all this groundwork, the student can now proceed to solve the TOV equations 
as before for a pure neutron star, using the fit for e{p) found in the previous sub-section. 
It is, once again, useful to convert from the units of MeV/fm^ to ergs/cm^ to Mo/km^ and 
dimensionless p and e. By now the student has undoubtedly grown quite accustomed to 
that procedure. 

e{p) = («:oeo)-'/y/' = Aop'/^ , Ao = 0.8642 , (98) 
where this time we defined 

4 5 

m„ c , . 

With this, the constant a that occurs on the right-hand side of the TOV equation, Eq. (j22I), 
is a = AqRq = 1.276 km. The constant for the mass equation, Eq. is P = 0.03265, 

again in units of 1/km^. 
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FIG. 12: The mass M and radius R for pure neutron stars using an EoS which contains nucleon- 
nucleon interactions. Only those stars to the right of the maximum are stable against gravitational 
collapse. Compare this graph with that in Fig. El which is based on a non-interacting Fermi gas 
model for the EoS. 

Now proceeding as before, one can solve the coupled TOV equations for p{r) and Ai{r) 
for various initial central pressures, p{0). We don't exhibit here plots of the solutions, as 
they look very similar to those for the Fermi gas EoS, Fig. El 

More interesting is to solve for a range of initial p(0)'s, generating, as before, a mass M 
versus radius R plot which now includes nucleon-nucleon interactions fFig. I12j). The effect 
of the nuclear potential is enormous, on comparing with the Fermi gas model predictions 
for M vs. R shown in Fig. IHl The maximum star mass this time is about 2.3 Mq, rather 
than 0.8 Mq. The radius for this maximum mass star is about 13.5 km, somewhat larger 
than the Fermi gas model radius of 11 km. The large value of maximum M is a reflection 
of the large value of nuclear (in) compressibility Kq = 363 MeV. The more incompressible 
something is, the more mass it can support. Had we fit to a smaller value of Kq we would 
have gotten a smaller maximum mass. 



E. What About a Cosmological Constant? 



We do not know (either) if there is one, but there are definite indications 
part of the make-up of our universe is something called "Dark Energy" 



that a great 
. This conclusion 
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comes about because we have recently learned that something, at the present time, is causing 
the universe to be accelerating, instead of slowing down (as would be expected after the Big 
Bang) . 

One way (of several) to interpret this dark energy is as Einstein's cosmological constant, 
which contributes a term Kg^^ to the right-hand side of Einstein's field equation, the basic 
equation of general relativity. The most natural value for A would be zero, but that may 
not be the way the world is. If A is non-zero, it is nonetheless surprisingly small. 

What would the effect of a non-zero cosmological constant be for the structure of a 
„eu.o„ .a., . tu.n. out ..at t.e ouly n,„diaca«o„ to t.e TOV equation Q . . the 
correction factor 



1 + 



M(r)c^ 



4'Kr^p{r) 



(100) 



M{ry 2GM{r)_ 

So, we encourage the student to, first, understand the units of A and then to see what values 
for it might affect the structure of a typical neutron star. 



VII. CONCLUSIONS 



The materials we have described in this paper would be quite suitable as an undergrad- 
uate thesis or special topics course accessible to a junior or senior physics major. It is a 
topic rich in the subjects the student will have covered in his or her courses, ranging from 
thermodynamics to quantum statistics to nuclear physics. 

The major emphasis in such a project is on constructing a (simple) equation of state. 
This is needed to be able to solve the non-linear structure equations. Solving those equations 
numerically, of course, develops the student's computational skills. Along the way, however, 
he or she will also learn some of the lore regarding degenerate stars, e.g., white dwarfs and 
neutron stars. And, in the latter case, the student will also come to appreciate the relative 
importance of special and general relativity. 
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[1] It is widely beleived that neutron stars were proposed by Lev Landau in 1932, very soon after 

the neutron was discovered (although we are not aware of any documented proof of this). 

In 1934 Fritz Zwicky and Walter Baade speculated that they might be formed in Type II 

supernova explosions, which is now generally accepted as true. 
[2] An on-line catalog of pulsars can be found at http:/pulsar. princeton.edu. 
[3] An intermediate-level on-line tutorial on the physics of pulsars can be found at 

|http : / / www. jb . man ■ ac. uk/ research / pulsar This tutorial follows the book by Andrew G. Lyne 
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lished in "The Nuclear Equation of State", ed. by A. Ausari and L. Satpathy, World Scientific 
Publishing Co., Singapore, 1996. 

[7] If you are a mentor for such a student, you may want to see some of the Mathematica and 
MathCad files we developed along our path of discovery. Send e-mail to the first author 
convincing him that you are such a mentor, and he will direct you to a web page from which 
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[9] S. L. Shapiro and S. A. Teukolsky, Black Holes, White Dwarfs and Neutron Stars: The Physics 
of Compact Objects, Wiley-Interscience, 1983. 
[10] We apologize to readers who are enthusiasts of SI units, but the first author was raised on 
CGS units. Actually, we strongly feel that by the time a physics student is at this level, he or 
she ought to be comfortable in switching from one system of units to another. 
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[11] A discussion of how to solve these equations (using conventional programming languages) is 
given in S. Koonin, Computational Physics, Benjamin-Cummings Publishing, 1986. 

[12] For more details on white dwarfs, NASA provides a useful web page at 
|http:/ /imagine.gsfc.nasa.gov/docs/science/know_ll / dwarfs.html 



[13] This maximum mass of 1.4 Mq is usually referred to as the Chandrasekhar limit. See S. Chan- 
drasekhar's 1983 Nobel Prize lecture, |http : / / www.nobel.se / physics /laureates / 1983 / For more 
detail see his treatise. An Introduction to the Study of Stellar Structure, Dover Publications, 
New York, 1939. 

[14] Mathematica is a software product of Wolfram Research, Inc., (see web page at http:/ /www.- 
wolfram.com), and its use is described by S. Wolfram in The Mathematica Book, Fourth Ed., 
Cambridge University Press, Cambridge, England, 1999. However, whenever we use the phrase 
"using Mathematica," we really mean using whatever package one has available or is familiar 
with, be it Maple, MathCad, or whatever. We did almost all of the numerical/symbolic work 
that we describe in this paper in Mathematica, but some of its notebooks were duplicated in 
MathCad, just to be sure it could be done there as well. 

[15] Enough of these explicit flags! Most of the equations from here on present challenges for the 
student to work through. 

[16] For the Newtonian case, a polytropic EoS also allows for a somewhat more analytic solution 
in terms of Lane-Emden functions. See Weinberg, op. cit.. Sec. 11.3, or C. Flynn, Lectures on 



Stellar Physics, especially lectures 4 and 5, at http : / / www. astro.utu.fi~cflynn / Stars /[ 
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[20] This fit is least accurate (~ 2%) at very low values of kp. However, this is where the pure 
neutron approximation itself is least accurate. The surface of a neutron star is likely made of 
elements like iron. A fictional account of what life might be like on such a surface can be found 
in Robert Forward's Dragon's Egg, first published in 1981 by Del Rey Publishing, republished 
in 2000. 
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[21] See Weinberg, op. cit., Sec. 11.2. 

[22] Because it is almost non-interacting with nuclear matter, a neutrino tends to escape from the 
neutron star. This is the major cooling mechanism as the neutron star is being formed in a 
supernova explosion. George Gamow named this the URCA process, after a Brazilian casino 
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studying reactions like Au nuclei striking each other at center of mass energies around 200 
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